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Abstract-In this article a new approach is considered for 
implementing multi-grid methods in iterative operator splitting 
methods for differential equations. The underlying idea is to 
embed fast multi-grid methods to accelerate the iterative splitting 
schemes. We concentrate on convection-diffusion-reaction 
equations and treat a multi-scale problem. The equations are 
spatially discretised with finite difference or finite volume 
methods. The main problem of iterative solvers is their neglect of 
multi-scale physics: different scales are not resolved on different 
spatial scales. Here we apply multi-grid methods to obtain the 
resolution of the optimal scale, which can be solved by an 
iterative splitting scheme. We discuss the embedding of such 
spatial- and time-scale methods and their underlying analysis. 
Numerical examples are given to verify the numerical schemes. 
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I. INTRODUCTION 

Our motivation is to solve multi-scale problems that occur, 
for example, in transport-reaction problems in porous media. 

In the last years, the interest in numerical simulations of 
multi-scale problems that can be used to model different 
physical behaviors, e.g., potential damage events, has 
significantly increased, but the adaptation of this methodology 
to such problems is delicate, see [39, 40]. 

We will concentrate on simplified models and give an 
account of the ideas needed to adapt the underlying splitting 
schemes, see [15]. While these splitting schemes worked well 
for single scale problems, we have to embed multi-grid ideas 
or multi-step ideas into these schemes to extend them to be 
multi-scale solvers. We will deal with iterative splitting 
methods and how to decompose such scale problems and 
accelerate the convergence rates with multi-scale methods, so 
that we can overcome such scaling problems. 

Moreover the embedding part is important: while 
dealing with multi-grid methods, we have to extend 
the underlying splitting analysis. We will discuss the 
needed algorithmic ideas. 

The novelty is to cheaply embed multi-grid and multi-step 
methods and apply only cheap iterative schemes to achieve 
higher order results. 

In general, we discuss an elegant way of embedding and 
survey recent methods for multi-scale problems with respect 
to iterative splitting schemes. 

In the following, we describe our model problem. The 
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model equation for the multi-scale equations are given as 
coupled partial differential equations. 

We concentrate on a far-field model for a plasma reactor, 
see [16] and [30], and assume a continuum flow and that the 
transport equations can be treated with a 
convection-diffusion-reaction equation, due to a constant 
velocity field: 

^ + V-Fc= 0, infix [0,t] (1) 
F = v - DV, 

c(x,t) = c (x),onfi, (2) 
0,ondftx [0,t] (3) 

where c is the particle density of the ionised species, F is 
the flux of the species, v is the flux velocity through the 
chamber and porous substrate, which is influenced by the 
electric field, and D is the diffusion matrix. The initial value is 
given as cO and we assume a Neumann boundary condition. 

The aim of this paper is to present a novel iterative 
splitting method that embeds multi-scale methods for spatially 
discretised transport-reaction equations. 

This paper is organized as follows: In Section II , we 
present the underlying splitting methods. The stability analysis 
of the embedded multi-grid schemes with the iterative 
schemes is given in Section III. The assembling and the 
algorithms for the embedded multi-grid method are given in 
Section IV. Numerical verifications are given in Section V. In 
Section VI, we briefly summarize our results. 

II. ITERATIVE SPUT™G„MErapDS 

Iterative splitting methods were developed in the early 90s, 
see [28, 41]. 

The idea is to decouple the equations into two or more 
equations to save computational time, without changing the 
previous results. We obtain inhomogeneous partial differential 
equations and solve them with the appropriate methods. 

We consider the following differential equation problem, 
where the operators A and B are spatially discretised: 

— = Au + Bu (4) 

where the initial conditions are un = u(tn). The 
solutions correspond via the spatial discretisation by u(t) 

= (ul (t), . . . , um (t))T = (c(xl , t), . . . , c(xm , t)T and m 
is the number of spatial discretisation points. The operators 
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A and B are matrices corresponding to the discretised 
convection and diffusion operators in (1) and their boundary 
conditions. Hence, they can be considered as bounded 
operators with a sufficient large spatial step Ax > 0. 

In the following we discuss the different methods based on 
the following definition: 

Definition 1. We define a sequential method as a serial 
method, while each step can only be done successively (e.g., 
each processor has to wait till the results of the previous 
computation are done). We define a parallel method as a 
multitasking method, which can be executed one piece at a 
time on many different processing devices, and then put back 
together again at the end to get the correct result, (e.g., each 
processor can independently compute a result). 

A. The Sequential Iterative Splitting Method 

The classical sequential operator splitting methods, known 
as Lie-splitting or Strang-Marchuk splitting, have several 
drawbacks besides their benefits, see [38, 31]. For instance, 
for non-commuting operators there might be a very large 
constant in the splitting error which requires the use of an 
unrealistically small time step. Also, splitting the original 
problem into different sub-problems with one operator, i.e., 
neglecting the other components, is physically questionable. 

In order to avoid these problems, one can use the iterative 
operator splitting method on the interval [0, T]. This algorithm 
is based on an iteration with a fixed splitting discretisation 

step-size x . On every time interval [t n , t n+ l] the method 
solves the following sub-problems consecutively for i = 0, 
2, . . . 2m. 



duj(t) 
9t 



= Aui(t) + Bui.^fhwith Ui (t n ) = u n ,t £ (t n ,t n+1 ), (5) 



= AuKt) + Bu i+1 (t),with u i+1 (t n ) = u n , t £ 
(t n ,t n+1 ), 



(6) 



where u n is the known split approximation at the time level t 

= t n (see [9]).. This algorithm is an iterative method which 
involves in each step both operators A and B. Hence, there is 
no real separation of the different physical processes in these 
equations. 

Remark 1. For the presented iterative splitting method, we 
have a serial algorithm and we cannot use parallel methods in 
this version. For the sake of efficiency, we will modify this 
method to enable parallelisation. 

B. The Parallel Iterative Splitting Method 

To parallelise the standard iterative splitting method, see 
Section II -A, we propose a decoupled version. 

The iterative scheme can be done in a Jacobian form, so 
that the two parts of the algorithm can be computed 
independently. 

In the next iteration steps, we use the previous solutions 
and improve them. 

We apply the parallel iterative splitting method for i = 0, 
2, . . . 2m, and we have: 



duj(t) 

at 



= Aui(t) + Bui.^t), with Uj(t n ) = u n , t £ (t n ,t n+1 ) (7) 
u (t n ) = u n ,u i _ 1 (t) = 0, 



3uj+i(t) 

at 



= AUj_ 2 (t) + Bu i+1 (t),with u i+1 (t n ) = u n ,t £ 



(t n ,t n+1 ), 
u_ 2 (t) = 



(8) 



where u n is the known split approximation at the time level 
t=tn. This algorithm is an iterative method where each step 
involves both Operators A and B. Hence, there is no real 
separation of the different physical processes. 
The first iterative steps are: 
au (t) 



at 

3ui(t) 

at 



= Au (t),with u (t n ) = u„, (9) 
= B Ul (t), (10) 



withUiCt") = u n , 



(11) 



and 



at 

du 3 (t) 



at 



= Au 2 (t) + BiiiOO.with u 2 (t n ) = u , 
= Au (t) + Bu 3 (t), 



with u 3 (t n ) = u n , 



(12) 



Remark 2. The effect of the parallelisation is to obtain the 
accuracy of the sequential iterative splitting method with one 
more iterative step. 

C. The Multi-grid Algorithm 

To describe the multi-grid method we first implement the 
two-grid method. Afterwards, we extend it recursively to the 
multi-grid method. The smoother grid Level 1 is denoted by 
SI. The two-grid method is defined as follows: 



Mf G := S^ 2 Mf GG S^ 1 



(13) 



where vl denotes the pre-smoothing steps and v2 denotes 
the post-smoothing steps. The correctionMf 00 on the coarse 
grid is defined by: 



Mf GG := 1 - pA^rAi 



(14) 



The passage to the multi-grid method is done in such a 
way that Al-1 of the coarse grid in Equation (14) is not 
inverted exactly, but the two-grid method is invoked y-times 
to solve the equation systems at the grid level 1—1. The 
system of equations is only solved on the coarsest grid. 



The multi-grid method is defined as: 



M™ u := 0, 
M, MG == M, ZG 



M mg 



Mf G 



+ S^ 2 p(M 1 N ! G ) y Ar_ 1 1 rA 1 S 



(15) 
(16) 
(17) 



where vl denotes the pre-smoothing steps and v2 the 
post-smoothing steps. The corrections Mf GG of the coarse 
grid are: 

Mf GG := 1 - p(/ - (MfLlY^yA^ (18) 

This approach is called the multi-grid method. For a 
choice of y =1 one speaks of a V-cycle, for y = 2 it is a 
W-cycle. 
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The multi-grid algorithm is given by 
Algorithm 1 Linear multi-grid cycle 

MG (x b b h l) 
{ 

m == 0) 
{ 

x = A^bf,; exact solving on coarse grid. 
} 

else 
I 

x 1 — S v± (xi, ftj); v 1 pre — smoothing steps. 
— rbi, defect restricted on 
next coarser grid. 

x h = 0; 

for G=0; i < y; 
{ 

ci-i = 0; 

MG(ci-\. bi-i, I - 1); y — fold recursive 

invoking of correction of coarse grid 
x l-i - x l-i + c i-l> 
Vi = bi-i - Ai_iCi_ 1; 

} 

q = pX[_!; interpolaton of correction 
of the next coarse grid. 

x, = X! + c,; 
b, = bi - Ajd; 

X[ = S V2 (x 1; b[); v 2 post — smoothing steps. 

} 

} 

The refinement is done using the strategy of [18]. By 
assuming a linear effort for the smoothers as well as for the 
grid-transfer operators, one obtains a linear effort for 
Algorithms 1, if y <3, cf. [18]. The proofs of convergence for 
the W-cycle were given in [17]. The topic of convergence will 
not be further discussed, for an overview we refer to [44]. 

Subsequently, the multi-grid cycles are illustrated in Fig. 1 . 
For a further treatment of the topic of multi-grid methods, we 
refer to [2, 17, 18]. 



Level 3 



■ Level 2 * 4 



■ Level 1 



J 

• * Voiglatteu 



Level 



A/W 



Nackglarten 

I Exakt Loseii 



V-Zyklus 



W-Zyklus 



Fig. 1 Multi-grid cycles: V- and W-cycle 

D. Iterative Splitting Method with Embedded Multi-grid 
Method 

The following algorithm is based on embedding the 
multi-grid method into the operator splitting method. The 
iteration with fixed splitting discretisation step- size t . On the 
time interval [tn , tn+1 ] we solve the following sub-problems 
consecutively for i = 0, 2, . . . 2m. (cf. [24] and [9] .) 



^ = AqW + P'a-'b + Bq-iCt), with Ci (t n ) = c n , (19) 

Sc i+1 (t) 



9t 



= R'A-'BAqCt) + Bc i+1 (t),with c i+1 (t n ) = c n , (20) 
where c (t n ) = c n , c—j = and c n is the known split 

approximation at t = t n . We assume A is the fine spatial 
discretised operator on level IA , where B is the coarse 
discretised operator on level IB . 

The operators are coupled by the restriction and 
prolongation operators: 



A — p'a-Ib a 

-"coarse 1 ' 
Bfine = P lA " lB B, 



(21) 

(22) 



where R is the restriction and P the prolongation operator. 

Algorithm 2 We assume two spatial operators A, B, which 
are discretised by finite difference of finite element methods. 



We solve the time-discretisation of our equations 



dq(t) 



9t 

^i±l™ = R'A-iBAqCt) + Bc i+1 (t), with c i+1 (t n ) = c n , (24) 



Aq(t) + P'A-'BBq.iCfhwith q(t n ) = c n , 



(23) 



where c0 (f 1 ) = c n , c—] = and c n and cn is the known 
split approximation at time-level t = f 1 by integration with 
respect to time and obtain 
(I - A)q(t n+1 ) = q(t n ) + P 1 A- 1 BBq_ 1 (t n+1 ), 

withq(t n ) = c n (25) 
(I - B)c i+1 (t) = q +1 (t n ) + R 1 * - ' 8 Aq (t n+1 ), 

withq +1 (t n ) = c n , (26) 

We have the multi-grid equations: 
Liq(t n+1 ) = q(t) n + P 1 A- 1 BBq_ 1 (t n+1 ), 

withq(t n ) = c n , (27) 
L i+1 q +1 (t) = q +1 (t") + R'A-'BAqCt^ 1 ), 

withq +1 (t n ) = c n , (28) 

E. The Waveform Relaxation Method 

A further method to solve large coupled differential 
equations is the waveform relaxation scheme. 

The iterative method was discussed in [41] and can be 
done either with Gaussian or Jacobian form. 

We deal with the following ordinary differential equation 
or assume a semi -discretised partial differential equation: 

ut =f(u, t), in (0,T), u(0) = vO 

where u = (ul , . . . , um)t and f(u, t) = (fl (u, t), . . . ,fm ( u > 

t)f. 

We apply the waveform relaxation method for i = 0, 1 , . . . m 
and have: 
du u (x,t) 



at 



du 2 ,i(x,t) 

at 



fl( u l,i- u 2,i-l- ■■■> u m,i-l) 
with U 14 (t n ) = Ui(t n ) 

= f2( u l,i-l< U 2,i' u 3,i-l' ■■■' u m,i-l) 
With u 2ii (t n ) =U 2 (t n ) 



(29) 



(30) 
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du mi (x,t) 

at 



= f 2 (ui,i- 



-l,i-l> u m,i) 



with u m ,i(t n ) = UnU (t n ) 



(31) 



where for the initialisation of the first step we put Uj _ y ( t) 
= u,(t n ),..., u^/t) = ujt n ). 

We reduce to two equations and reformulate the method to 
our iterative splitting methods. 



So we consider 



at 

du 2 



= f 11 (u 1 ,t) + f 12 (u 2 ,t),in (0,T) 
(t =f 21 (u 1 ,t)+f 22 (u 2 ,t),in(0,T) 
u(0) = v 

where u = (uj ,U2 f. 
Our notation for the operator equation is 
g = A(u) + B(u),in (0,T) 
u(0) = v 

Where 



(32) 

(33) 
(34) 



B(u) = (["^ 



(Ml) J 
(u 2 )\ 

(Hi)) 



(35) 
(36) 

(37) 
(38) 



The iterative splitting method as waveform relaxation 
method written for 

i = 0, 1, . . . m as: 

-^i = A(u li ,u 2i _ 1 ) + B^j.i.Uy) 

with u^Ct") = Ui(t n ) and u 2 i (t n ) = u 2 (t n ) 



(39) 



where for the initialisation of the first step we put m ; _ 1 (t) = 

U l U 2,-l(t) = U 2 (t U ). 

Convergence analysis In the following we formulate the 
iterative splitting method as a waveform relaxation method 
and apply the proof techniques of those relaxation schemes. 

The waveform relaxation method is given by 



diij 
dt 



P Ui + QUi.! + f, 
Ui (t n ) = u(t n ), 



(40) 
(41) 



where A = P + Q, e.g., P is the diagonal part of A (Jacobi 
method). 

Here the splitting method is done abstractly with respect to 
the Matrix A. The method is considered an effective solver 
method with respect to the underlying matrices. 

The reformulation of the iterative operator splitting 
method is given by 



diij 
dt 



= P Ui + QUi.! + f, 



Ui (t n ) = u(t"), 



dui. 



dt 



P Ui + Qu i+1 + f, 



(42) 
(43) 
(44) 



u i+1 (t") = u(t n ), 



(45) 



where P, Q are matrices given by spatial discretisation, e.g., 
P is the convection part of Q the diffusion part. 

But we can also make an abstract decomposition, taking 
into account that A = P + Q, where P is the matrix with small 
eigenvalues and Q is the matrix with large eigenvalues. 

Theorem 1 We have the iterative operator splitting 
methods given as the following outer and inner iterations. 



Outer Iteration: 



dUj 
~dT 



= PUi + QUi.! + F, 



Ui(t n ) = U(t n ), 



(46) 
(47) 



where P and Q are the diagonal and outer-diagonal 
matrices of the splitting methods. 



Inner Iteration: 



du 
dt 



I _ 



P Uj + QUj.! + f, 

Uj (t n ) = u(t n ), j = l J 

^ = Pu, + Qu k + f, 
u k (t n ) = u(t n ), k=l K 



(48) 
(49) 
(50) 
(51) 



We have a convergent scheme for K bounded in a 
Banach-space: 



p(K) : = lim ||K K || 1/K 



(52) 



where p(K) < 1 is the spectral radius of K and is given by 
the variation of constants of P and Q. 



Proof. The outer iteration is given by 

^ = PUi + QU^ + F, 
Ui(t n ) = U(t n ), 



(53) 
(54) 



where P and Q are the diagonal and outer-diagonal 
matrices of the splitting methods, which can be solved by the 
variation of constants. 



We introduce the linear integral operator 
KU(t) = J t exp((t-s)P)QU(s)ds, 



And we have 

U, = KUi.! + 0, 



(55) 



(56) 



= exp(tP)U(t n ) + J t exp((t - s)P) F(s)ds, (57) 

For the convergence it is sufficient to show that K is 
bounded in a Banach-space: 

p(K) := lim ||K K || 1/K (58) 

K— >co 

where p(K) < 1 is the spectral radius of K 

This is given by the bounded operators, see [33, 34]. 
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III. ERROR ANALYSIS: STABILITY THEORY FOR THE ITERATIVE 
SPLITTING METHOD WITH MULTI-GRID METHODS 

The following algorithm is based on embedding the 
multi-grid method in the operator splitting method. The 
iteration with fixed splitting discretisation is step-size x. On 

the time interval [t n , t n + 

sub-problems consecutively for i = 0, 2, . . . 2m. (cf. [24, 9].) 

^ = Aq(t) + P'A-'BBq.iCtXwith Ci (t n ) = c n , (59) 
^i±i™ = rU-Ibac.^) + Bc i+1 (t),with c i+1 (t n ) = c n , (60) 

where c g (t n ) = c n , c_ ; = and c n is the known split 

approximation at time-level t = f 1 . We assume A to be the 
fine spatial discretised operator on level IA and B to be the 
coarse discretised operator on level IB. 

The operators are coupled by the restriction and 
prolongation operators: 



A - [r1a-1b°a b ] 



(66) 



Then using the Notation (66), the Relations (64) and (65) can 
be written in the form 



we solve the following a t£i (t) = A £i (t) + F^tXt £ (t n ,t n+1 ],£i(t n ) = 0. (67) 



Due to our assumptions, A is a generator of the one-parameter 
Co-semi-group (exp At) t > 0, hence using the variations of 
constants formula, the solution to the abstract Cauchy problem 
(67) with homogeneous initial condition can be written as 

£i (t) = f exp(A(t - s)) Fi(s)ds, t E [t n , t n+1 ] (68) 



A — rU-Iba 

^coarse Iv "i 
Bfine = P lA " lB B, 



(61) 
(62) 



where R is the restriction and P the prolongation operator. 

Theorem 2 Let us consider the abstract Cauchy problem 
in a Banach space X 



d t c(t) = Ac(t) + P'a-^BcOO, < t < T 
c(0) = c o 



(63) 



where A, P ,A ~ lB B, A + P lA ~ lB B : X 



X are given 

linear operators being generators of the C -semigroup and c 
£ X is a given element. Then the iteration process (59)-(60) 
is convergent and the rate of convergence is of higher order. 

Proof. We assume A + P ,A ~' B B <E L(X) and assume a 
generator of a uniformly continuous semi -group, hence the 
problem (63) has a unique solution c(t)=exp((A + P 
B)t)c . 



I A -IB 



Let us consider the Iteration (59)-(60) on the sub-interval 
[tn, tn+1]. For the local error function ei(t) = c(t) - ci (t) we 
have the relations 

a t ei(t) = Aei(t) + P'A-'BBei.^t), t £ (t n ,t n+1 ] 

ei (t n ) = 0, (64) 

And 

d t e i+1 (t) = R'A-'BAeiCt) + Be i+1 (t), t £ (t n ,t n+1 ] 

e i+ i(t n ) = 0, (65) 

for m =0, 2, 4, . . ., with e0(0) = and e~l(t) = c(t). In the 
following we use the notations X2 for the product space X x 
X endowed with the norm 

||(u,v)|| = max{||u||, ||v||}(u,v £ X) The elements 

£i(t). F[(t) £ X 2 and the linear operator 

A : X 2 — > X 2 are defined as follows: 



ei(t) 
e i+ i(t) 



,Fi(t) = 



piA-iBBe^^t) 




(See, e.g., [8]) Hence, using the notation 

llEjlloo = SUp te [ t n t n + l]||£i(t)|| (69) 

We have 

l|Eill(t)< ||F i |L/ t t n ||exp(A(t-s))||ds = 
||B||||e i _ 1 ||/ t t n ||exp(A(t-s))||ds,t£ [t",t n+1 ]. (70) 



Since (A(t))t > is a semi-group, the so called growth 
estimation 

1 1 exp (At) 1 1 < Kexp(wt),t > 0, (71) 



holds with some numbers K > and co £ IR, cf. [8]. 

Assume that (A(t))t > is a bounded or an exponentially 
stable semi-group, i.e., (71) holds with some co < 0. Then 
obviously the estimate 



I exp (At) 1 1 < K,t > 0, 



(72) 



holds, and hence, on the basis of (70), we have the relation 

IN|(t) < K||P 1 A-iB B ||T n ||e I _ 1 || > tE [t",t n+1 ] (73) 

- Assume that (exp At)t > has exponential growth with 
some co > 0. Using (71), we have 

f||exp(A(t-s))||ds<K w (t),t£ [t",t n+1 ] (74) 



Where 



K w (t) = -(exp(w(t - t n )) - l),t £ [t n ,t n+1 ] (75) 
w 



Hence 



K w (t) < -(exp(wT n ) - 1) = K Tn + 0(T 2 ). 



(76) 



The Estimations (73) and (76) result in 

llqlloc = KlIP'A-^HllBIITJIe^ll + 0(T 2 ) (77) 

Taking into account the definition of Ei and of the norm k ■ koo, 
we obtain 

||e,|| = KlIP'A-'BllHBllTJIei.J + 0(T 2 ) (78) 
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And hence 

lle, +1 || = K 1 T n 2 ||e 1 _ 1 || + OCr„ 3 ). 
which proves our statement. 



(79) 



Operator-splitting method with embedded Jacobian 
Newton iterative method: Newton's method is used to solve 
the nonlinear parts of the iterative operator-splitting method, 
see the linearisation techniques in [24, 25]. We apply the 
iterative operator-splitting method and obtain 

FiCuj) = d t Uj - A(Ui)Ui - BO^Ou^ = 0, 

with Ui (t n ) = c n , 

FzOi+i) = d t u i+1 - A^H - B(u i+1 )u i+1 = 0, 
withu i+1 (t n ) = c n , 

where the time step isr = t n + ^ — t n . The iterations 
are i = 1, 3, ... , 2m + 1. c (t) = is the starting solution and 
cn is the known split approximation at time level t = t". The 

results of the methods are c(t n + ^ ) = u2m + 2 (t n + ^. The 
splitting method with embedded Newton's method is given by 

u[ k+1) = u[ k) - D (F^uf ))" 1 (a t u[ k) - A( Ul (k) )u[ k) 
-B(u I (k) 1 )u I (k) 1 ),withD(F 1 (u I (k) )) 

._( A(u « )+ ^p u «), 

and k = 0,1,2, ... , K, with Uj(t n ) = c n , 

-l 



(k) 



u^=u»-D(F 2 (u«))- (a t u»-A(u[>[ 
-B(u[ k l)u«)c»), 

with D (F 2 (u»)) = - (b(u«) + ^f^uffA 



du\ +1 

andl = 0,1,2, L, with u ; + l(t n ) = c n . 

Remark 3. For the iterative operator-splitting method with 
Newton's method we have two iteration procedures. The first 
iteration is Newton's method for computing the solution of the 
nonlinear equations; the second iteration is the iterative 
splitting method, which computes the resulting solution of the 
coupled equation systems. The embedded method is used for 
strong nonlinearities. 

IV. ASSEMBLING OF THE SPLITTING METHODS WITH THE 
EMBEDDED MULTI-GRID METHOD 

In the following we discuss the assembling of our methods. 

A. Crank-Nicolson and Two-grid Method 

In the following we apply the iterative splitting method, 
with a Crank-Nicolson method in time, a second order finite 
difference method in space, and an underlying two-grid 
method as solver. 

We apply the Crank-Nicolson scheme separately to the 
two equations obtained by the two steps of the iterative 
operator splitting method after the discretisation of the 2D 
heat equation. For the ith step of the iterative operator splitting 
method, as interpolation operator we choose the bilinear 
interpolation given by 

h _ 2h 
v 2i,2j — v ij 



v 2 h i+1 , 2j = 2K h + v ^.i)' 

V 2 h i,2j + l=^(vf +V?f +1 ) 
1 

v h _ t ( 2h , 2h , 2h , 2h ^ 

v 2l+l,2j + l — ^ ^ V 1J T v 1+1) T V 1J + 1 T V l+1] + 1 J, 

0<i,j,<^-l. 

where the superscript h is associated with the fine grid and 2h 
with the coarse grid. 



Step i: 



-H.j.k _ u i,j,k 



At 



u i-l,j,k zu i,j,k + u i+l,j,k 



h 2 



+ PD, 



u i-l,j-l,k-l zu i-l,j,k-l ^ u i-l,j+l,k-l 



h 2 



+ PD 2 



u P-l,j,k _ 2u"j k + uf +1 j k 



h 2 



<^> Puj 



n+1 



U-^-U 2h 2 



u f-ij-i,k-i _ ^ u f-i,jM-i + ujLij+^k-iN 
h| ) 

,j,k-l + 2 °x u i,j,k + 

-D 2 At\ n+1 /-DjAt 



D 7 At 



n+1 

,j,k 



+ <k^ + Pu n+1 



U,k h 2 



i-l,j,k-l . 2 

Hy 

-D,At 



2h 2 



+ 



u n + l (Zh^\ 
Ui+1 H 2h? ) 



= 2 ( D i 5 x<J,k + D 2 P5 y ur_ 1 , j , k _ 1 ) + < k 



<=> u 



n+1 / -DzAtA 

2h y ) + u ^H 2h 2 J + Ui 'i 

+ u n + l^ + ir u n + l +u n + l ^2* 
+ U U.k ^2 + 2 V U i-l,j-l,k + u i-lj+l,kJ 



+1 
j,k 



h 2 



+ ^ { U i-1,)+1M + U U+l,k + u Kj+2,k + u Sj+2,k) y 



+ 



u n + l 

Ui+1 H 2h 2 ) 



-D 2 At 



= -D 1 5X jk + -D 2 P 
+ u "j,k 



1 u i-l,j-l,k-l U i-1 ,j+l,k-l 



+ ^ ( u ?-l,}+l,k + u ij+l,k + u ?-l,]+2,k + u Sj+2,k) ^ 
+ u n + l 

+ Ui+1 'i' k V 2h 2 J 

_i 2n id 2 r n i 

= 2~Di8 x u]j k + 2 h 2 ^ l^-y-Lk _ ^ - (ujLij-^k + up_ lj+lk ) 



D 2 At 

-D 2 At 
2h|" 



+ ■ 



u 



i-l,j+l,k 



4 \+ u S!j+l,k + u P-l,j+2,k + u i!j+2,k 

n /-DiAt\ _,,/ Dt At 



/D 2 At> 
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+u 



-D,At 



-D,At 



Step i + 1: 



D 2 At 



+<-ti+2,k^ 4h 2 



+ 



u n + l f -^ At ) 

Ui+1 -i- k V 2h? J 



ll n+1 _ ,, n 

u i+l,j,k+l u i+l,j,k+l 



At 

(n n+1 _ 7i, n+1 4- n n+1 
u i-l,j,k zu i,j,k + u i+l,j,k 
RDl hi 

,,n+l n,,n+l I ,,n+l 



+ T ( u i!j+l,k + u i'-l,j+2,k + u i!j+2,k) 



DtAt D 2 At 
h\ ' Ihl 




f \ 
-D 2 At 
ihl 


■ „n+l 
+ U i-1 ,j,k 


I ) 








( 


\ 




\ 



+ u u,k 

DiAt 

2h\ 



j_>j n+1 



-D 7 At 



+ u s+i,i,k 



-DjAt 
1 



i, n+1 _ 7ii n+1 4- n n+1 \ 

n u i+l,j-l,k+l iu i+l,j,k+l T u i+l,j+l,k+l \ 

/ u i-l,j,k ~~ 2u i j k + U i+l j k 

' n \ „ 

_ u i+l,j-l,k+l — 2u i+1 j k+1 + U i+1 j +lk+1 \ 

+ D2 ^ ) 

— RK 2 ,, n+1 -I- ^ 2 K 2 ,,n+1 _i_ RX2,,n 

- ~^~ K °x u i,j,k + ~^~°y u i+l,j,k+l + ~^~ K °x u i,j,k 

+ ~^^Y U ?+1 ,j,k+l 



<=> auij_ 2>k + bui_ lij>k + cu 1Jjk + bu 1+1J(k + au iJ+2jk = d n jk 
This procedure leads to a linear system Au = f , where 

u = [ u 2,2 u 3,2 ■ ■ ■ u n/2-l,2 > U 2,3 U 3,3 ■ ■ ■ u n/2-l,3 > ■ ■ ■ 
U2,n/2-lU3,n/2-l • ■ • Un/2-l,n/2-l ]T 

f = [d2,2 ^3,2 ■ ■ ■ d n /2-l,2 > d 2j 3 d 3 3 . . . dn/2-1,3 , ■ ■ ■ d 2i n/2-l 
d3,n/2-l • • • dn/2-l, n /2-l ]T 

The index k + 1 is omitted in the vectors u,/for the sake of 
better presentation. The matrix A is given by 



cb Q 0000 a 
b c b a 
Ob c b a 
b c b 
b cb 
be b 

b c b 

a b c b 

Da b c b 



b c b 



b c b 

a00000&c 



For the (i + l)st step of the iterative operator splitting 
method, as restriction operator we choose the most obvious 
restriction operator, the injection. We prefer it instead of the 
standard full weighting operator because it leads to a linear 
system with better structure that is easier to handle. Injection 
is defined by 

v 2h _ 2h 
v i,j — v 2i,2j- 



z 

n+1 /-°2 At \ n+1 /- D iAt\ n+1 

<+l,j-l,k+l \^r) + RU ""Uk (-^2-) + U P+l,j,k + l 

^ D x At A , D,At 



+ KU i.j.k h 2 + U i+l,j,k+l h 2 

/-D 2 At 

+ u n+l / 2 



■M+lj+lJfc+l 



2h 2 



j /-D.Atx 



= -(D 1 R5lu n jik + D 2 5 2 < +1 ,j, k+1 ) 



n n+1 

^ u i+l,j-l,k+l 



+ u f+l,j,k+l 

-D 2 At 



...n + 1 / -DlAt \ 

2h 2 ; + u i-Uk+iV 2h 2 J + Ui+ u 

+ U i,j,k+1 h? 



,j,k+l 



■ + uf 



4. ,,n+i 

+ U i+ l,j + l,k+l \ 2h 2 



l i+l,j,k+l ' 

-D 2 At> 



h 2 



4- li n+1 

T u i+l,j,k+l 



2hJ 



= ^(D 1 6 2 < jik+1 + D 2 5X+ij,k+i) 



u i+l,j,k+l 



** u i.j-l.i+l 



, ,.B+1 







( \ 






-D\At 
■Ihi 














f 


V 



„ J . 1 D|4* 



2^ 



— / 



.t+i 



<=> au u _ ljk+1 + bUi_ 1Jjk+1 + 
+ au 



cu 



i,j + l,k+l — d"j k+ i 



i,j,k+l + bu i+1 j k+1 



Au = f , where 



Un/2-1,3 , 



U = [U 2 ,2 U 3 , 2 . . . Una-1,2 , U 2>3 U33 
u 3,n/2-l ■ ■ ■ u n/2-l,n/2-l ]T 

f = [d2,2 &3,2 ■ ■ ■ dn/2-1,2 > d 2 ,3 d 3 3 . . . dn/2-1,3 ) 
d3,n/2-l ■ ■ ■ dn/2-l,n/2-l ]T 



U2,n/2-l 
d2,n/2-l 



(The index k+l is omitted in the vectors u, /for the sake 
of better presentation.) 
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.4 = 



cbd a 
b •■ b a 

6c b 

b c 

a b 
a 



a 

b 

c b 
b c b 
b c 



a b c b 

a b c 



The restriction operators are given by 

RiC 1 lead to i K_tj-i,k-i + 2u$_\ M _ x + uf+y-yc-i + 

n,,n+l _i_ /i,,n+l 4- ?n n+1 4- n n+1 _i_ 
zu i-l,j,k-l ^ u i,j,k-l zu i+l,j,k-l + u i-l,j+l,k-l 

? n+1 , n+l ] 
^ u i,j+l,k-l T u i+l,j+l,k-lJ' 

according to the interpretation of the stencil 



1 

16 



12 1 
2 4 2 
12 1 



V. NUMERICAL EXPERIMENTS 

In the following examples, we deal with different test 
examples and their underlying multi-scale physics. 



The Jacobi iterations for the heat equation are given by 



u[" +1) =AtD r 



u i-l,j T iu ij 



,(n+l) 
i+ij 



+ AtDz +Zu * ^±1. + u w 



h 2 



i] 



B. Implementation of the Iterative Splitting Method with 
Embedded Two -grid Method 

The implementation for the heat equation is given by 



Mj,k 



^i.j.k 



At 



u 



n+l 
i-l,j,k 



n n+l , n+l 
zu i,j,k + u i+l,j,k 



+ PD, 



u 



n+l 

i,j-l,k-l 



- 2u 



n+l 
i,j,k-l 



+ U 



n+l 

i,j+l,k-l 



+ D,R- 



u 



i-l,j,k 



h 2 

Uy 

2 u i,j,k + u i+l,j,k 



+ D ; 



u i,j-l,k+l ^ u i,j,k+l T u i,j+l,k+] 

h 2 

Uy 



- ^°x u i,j,k + ^™y u ij,k-l + ^T Ko x u i,j,k 



+ — S 2 n n 



<=> p<-i,k-l 



-D 2 At 
2h 2 



-D, At 



2h 2 



) + <$ 



+ u 



D,At 



i.j.k ~~u2~ + Pu U.k-l 

-D 2 At 



D 2 At 
h 2 

Uy 



I p,,n+l 
+ ru i,j+l,k-l 



2h 2 



/-D.,At\ 



+ u 



i,j,k 



(D^S^u^ + D 2 P6 y < jik+1 ) + < jjk 



\ 



+ K + hi ' 



7 



0" 



«* a <-i,k + buP+Vi + cu^ 1 + au^k + buP + tj, k +i 



(1) 



A. Simple Heat Equation 

We deal with a PDE which is parabolic and has stiffness in 
the diffusion part. 

We have the following equation: 

3 t Ui = D n + D 21 d xx u 2 , 

d t u 2 = D 12 d xx u ± + D 22 d xx u 2 , (80) 

u ± (0) = u 10 ,u 2 (0) = u 20 (initial conditions), (81) 

u x (x,t) = u 2 (x, t) = (boundary conditions), (82) 

where Dll , D21 , D12 , D22 E IR+ and Dll , D12 < 
D21, D22 are the diffusion operators. 

We apply the finite difference scheme for the spatial 
derivatives: 

5xx = [1-21] and 5yy = [1 - 2 l]t . 

With the standard projection and restriction: 

1 2 11 
R H = 1/16 2 4 2 
.1 2 1. 

Ph = 4R H 

The time-derivations are discretised by implicit Euler 
methods and we obtain the following system of linear 
equations: 

( Cl n+1 - cf)/At = Aq n+1 + PBcBSwith q(t n ) = c n (83) 



(c^ 1 - cf +1 )/At = RA Cl n+1 + Bcgi 1 . 

with q + l(t n ) 



(84) 



We obtain, after the insertion of the operators, a system of 
linear equations: 



AUi n+1 = BU^ 1 + CU? 1 , 



(85) 



This is solved by a two-grid method. Here we have two 
scales and decouple: 



d t u = Au + Bu, 
u(0) = (u 10 ,u 20 ) T , 

where u(t) = (u] (t), u2 (t)) T 'for t E [0, T]. 
Our split operators are 



(86) 
(87) 
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/D, 



^ = . 'ii 3 XX 0\ /0 D 21 d xx \ 

i 12 a xx oJ' Vod 22 5 xx , 



(88) 



We chose such an example to have AB #= BA, and 
therefore we have a splitting error of first order for the usual 
sequential splitting methods, called A-B splitting. 

Our numerical results based on two-grid methods are 
presented in the following Table I . 

We choose Dl 1 = D12 = 0.5 and D21 = D22 = 0.05 on the 
time interval [0, 1]. As second order method we choose 
Crank-Nicolson with 9 = 1/2. As fourth order method we 
choose the Gauss Runge-Kutta. We apply a two-grid method, 
see subsection IV-A. 

The numerical results are given in Table I . We apply the 
L2 error based on the exact and the numerical solution. 

TABLE I NUMERICAL RESULTS FOR THE FIRST EXAMPLE WITH THE ITERATIVE 
SPLITTING METHOD AND SECOND AND FOURTH ORDER METHODS 



more optimal results for the iterative schemes. 
The result is given in Figs. 2-4. 




30 35 




Fig. 2 Blue - Iterative Operator Splitting (4 iterations); Green - A-B Splitting 



Number 


Iterative 


em 


ETT2 


ern 




of time- 


Steps 


(order 2) 


(order 2) 


(order 4) 


(order 4) 


paxtitions 












2 


1 


4.5321e-002 


3.6077e-OC3 


4.5321e-O02 


3.6077e-003 


2 


10 


3.9664e-003 


4.7396e_004 


3.9664e-O03 


4.7397e-004 


2 


100 


3.92O4e-O04 


4.8078e-005 


3.9204e-0(M 


4.8C83e-005 


3 


1 


4.6l26e-004 


3.6077e-003 


4.S12&6-004 


3.6077e-003 


3 


10 


7.8129e-006 


2.9285e-005 


7.S0S9e-006 


2.9289e-005 


3 


100 


8.5988e-008 


2.8270e-007 


8.0050e-008 


2.8682e-007 


4 


1 


4.S12&6-004 


2.2459e-0C5 


4.6l2&e-O04 


2.2464e-005 


4 


10 


4.1S83e-007 


4.2629e_008 


4.1321e-O07 


4.8154e-008 


4 


100 


5.952le-009 


5.4846e-C09 


4.0839e-OlO 


4.9968e-0ll 



Remark 4. In the experiments, we improved the iterative 
splitting scheme with higher order discretisation methods. 
Further we can accelerate the splitting scheme with an 
underlying two-grid method that is embedded in the splitting 
scheme. Both higher order time-discretisation and also 
two-grid methods are necessary to obtain such results. 

B. Second Example: Transport— Reaction Equation 
We deal with the following system of transport equations: 

d t u ± + v ± d x u ± — — + uu 2 , (89) 

a t u 2 + v 2 d x u 2 = uui - }u 2 , (90) 



where we have Dirichlet boundary conditions and wi o • 
are the initial conditions. 



«2,0 



We split the operator into two operators: 



(91) 



We see that for (j. ~ the operators commute. 

We choose vl = 1, v2 = 0.5, X = 0.01. We let t e [0, 40] and x 

e [0,40]. 

We set 

1 4 
At = — and Ax = — 
25 10 

In Fig. 2, we choose u = 0.001 and see that we could say u 
~ and obtain nearly the same results as for the A-B splitting. 

In Fig. 3, we choose u = 0.01 and see that u = and obtain 




s m 




1= 3D 2 


5 4 




\ \ 







Fig. 3 T=2: Blue - Iterative Operator Splitting (4 iterations); Green - A-B 
Splitting 




Fig. 4 T=60: Blue - Iterative Operator Splitting (4 iterations); Green - A-B 
Splitting 

Remark 5. In this example, we concentrate on the 
comparison between the standard A-B splitting and the 
iterative splitting methods. We obtain more accurate results 
for non-commuting problems. Such problems are related to 
multi-scale problems and we can apply our embedded 
multi-grid method. We have the benefit of receiving higher 
order results without reducing the time-step. 

VI. CONCLUSIONS AND DISCUSSION 

We presented iterative splitting methods with embedded 
multi-grid and multi-stepping methods. Here the idea was to 
achieve more accurate and faster results than the standard 
splitting schemes (e.g., the A-B splitting schemes). We 
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obtained a benefit from this embedded method in accelerating 
the iterative schemes and achieved much more accurate results. 
In the future we will concentrate on nonlinear and matrix 
dependent splitting schemes. 
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